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ABSTRACT 

We present a new code designed to solve the equations of classical ideal magneto- 
hydrodynamics (MHD) in three dimensions, submitted to a constant gravitational 
field. The purpose of the code centers on the analysis of solar phenomena within the 
photosphere-corona region. We present ID and 2D standard tests to demonstrate the 
quality of the numerical results obtained with our code. As solar tests we present the 
transverse oscillations of Alfvenic pulses in coronal loops using a 2.5 D model, and as 
3D tests we present the propagation of impulsively generated MHD-gravity waves and 
vortices in the solar atmosphere. The code is based on high-resolution shock-capturing 
methods, uses the HLLE flux formula combined with Minmod, MC and WEN05 re¬ 
constructors. The divergence free magnetic field constraint is controlled using the Elux 
Constrained Transport method. 
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1 INTRODUCTION 

In solar physics one of the most important subjects is re¬ 
lated to the dynamics of the plasma in the solar atmo¬ 
sphere. Specifically, phenomena which are related to the 
propagation of magnetohydrodynamic waves (e. g. Alfven, 
slow, fast magnetoacoustic), and transients (e. g. vortices, 
jets, spicules). Recent observations based on data provided 
by current space missions, such as Transition Region and 
Coronal Explorer (TRACE), Solar Dynamics Observatory 
(SDO), Solar Optical Telescope (SOT), X-ray Telescope 
(XRT) and Swedish Solar Telescope (SST) ([Banerjee et al. 


2007| [Tbmczyk et al. 2007| |De Pontieu et al. 2007| |Okamoto 


et al. 2007| [Jess et al. 2009 ) have enhanced the study of 
these phenomena in the solar atmosphere. 

The type of effects in the solar atmosphere are con¬ 
sidered when the plasma is fully ionized and therefore the 
ideal MHD equations represent a good model. However in 
this limit, these equations are highly non-linear and numeri¬ 
cal simulations are essential, where accurate numerical algo¬ 
rithms play an important role. Over the recent years a num¬ 
ber of codes have been used to understand various effects in 
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the solar atmosphere. For instance, in (Konkol & Murawski 


2010 Danilko et al. 2012) the authors used a modification 


of ATHENA code (Stone et al. 2008) to study vertical os¬ 


cillations of solar coronal loops and magnetoacoustic waves 
in a corollarygravitationally stratified solar atmosphere. In 


( Konkol et al. 20Tol Murawski & Musielak 2010 Murawski 


et al. 2011 Murawski et al. 2013| Chmielewski et al. 2014a) 


the authors use a modification of the FLASH code (Fryxell 


et al. 2000 Lee & Deane 2009 Lee 2013) to simulate the 


propagation of MHD waves in two and three dimensions. 


Moreover in (Woloszkiewicz et al. 2014) the authors study 


the propagation of Alfven waves in the solar atmosphere 


using the PLUTO code (Mignone et al. 2007). 


Wave phenomena in the solar atmosphere have influ¬ 
enced the development of numerical codes dedicated specif¬ 
ically to study MHD wave propagation, e.g. the VAC code 


(Toth 2000 Shelyag et al. 2008), which is capable to sim¬ 


ulate the interaction of any arbitrary perturbation with a 
magnetohydrostatic equilibrium background. The SURYA 


code (Fuchs et al. 2011), allows to study the propagation 


of waves in a stratified non-isothermal magnetic solar atmo¬ 
sphere, in particular in the chromosphere-corona interface. 
Numerical simulations also permit to analyze some transient 


effects, for example in (Murawski & Zaqarashvili 2010 Mu¬ 


rawski et al. 2011) using FLASH, the authors simulate the 


formation of spicules and macrospicules in the solar atmo- 
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sphere. In addition in (Pariat et al. 2009 Pariat et al. 2010 


Pariat et al. 2015) the authors propose a model for jet ac¬ 


tivity that is commonly observed in solar corona. Also the 
3D ideal MHD equations are used to model the magnetic re¬ 


connection producing massive and high-speed jets (DeVore 


1991). 


Even though all the aforementioned codes are available, 
well documented and can be applied to various problems, 
we present our new application, with convergence tests and 
error analyses. Our code is at the moment useful to study 
solar phenomena in the photosphere-corona region, however 
with plan to generalize it to deal with magnetic reconnection 
and solar jets including radiative processes. 

The paper is organized as follows. In Section we de¬ 
scribe the ideal MHD equations and mention the numerical 
methods implemented in our code. In Section we show 
the standard tests in ID and 2D dimensions. As a solar 
test in 2 D we study the transverse oscillations of Alfvenic 
pulses in the solar coronal loops. We also present as a 3D 
test the propagation of MHD waves in a gravitationally non- 
isothermal solar atmosphere. Finally, in Section^ we men¬ 
tion some final comments. 


magnitude g — 274 m/s^, which is the value of the con¬ 
stant gravitational acceleration on the solar surface. These 
equations are written in units such that /io = 1- It is also 
convenient to write the equation of state in terms of the 
temperature, that is 


rr 


( 8 ) 


where ks is Boltzmann’s constant, m is the mean particle 
mass defined in terms of the mean atomic weight g in the 
following way g = mimp, where nip is the proton mass. The 
above mentioned system of equations 0-0 can be written 
as a first-order hyperbolic system of flux-conservative equa¬ 
tions, which in cartesian coordinates reads 


dQ dF dG du 

dt dx dy dz 

V - B = 0, 

with Q a vector of conservative quantities 


(9) 


2 MHD EQUATIONS AND NUMERICAL 
METHODS 


Q = {p^pVx,pVy,pv^,E,B:,,By,B^)^ , 


2.1 Ideal MHD Equations 

We consider the set of the ideal MHD equations submitted 
to an external constant gravitational field. We choose to 
write down these equations in a conservative form, which is 
optimal for our numerical methods: 


^+V.(pv) = 0, 

+ V-(pv(g)v-B(g)B + Ipt) = pg, 

-^ + V • ((S + Pt)v - B(v • B)) = pv • g, 

dB 

^ + V-(v(g)B-B(8)v) = 0, 

V - B = 0, 


( 1 ) 

( 2 ) 

(3) 

(4) 

(5) 


where p is the plasma mass density, pt — P + ^ is the to¬ 
tal (thermal + magnetic) pressure, v represents the plasma 
velocity field, B is the magnetic field, I is the unit matrix 
and E is the total energy density, that is, the sum of the 
internal, kinetic, and magnetic energy density 


and F, G and H are the fluxes along each spatial direction 


F = 


pVx 

pvl - Bl+pt 

pVxFy BxBy 
pVxVz - BxBz 

{E-\-pt)vx - Bx{B-v) ’ 
0 

^xBy VyBx 
VxBz - VzBx 


G = 


pVy 

PVyVx ByBx 

pvl - Bl+pt 

pVyVz ByBz 
{E E Pt)vy - ByiB • w) ’ 
^yBx XxBy 

0 

'^yBz VzBy 


E = 



(6) 


The internal energy term is calculated considering an equa¬ 
tion of state of an ideal gas 


P = (7 - l)pe, (7) 

where e is the internal energy density and 7 is the adia¬ 
batic index. This equation also serves to close the system 
of equations 0 - 0 . In solar surface problems one considers 
the gravitational source term given by g = [0, 0, —g] with 


H = 


pVz 

pVzVx - BzBx 

pVzVy - BzBy 

pvl - Bl+pt 
{E-\-pt)vz - ^^(B • v) 

'IdzBx VxBz 
VzBy - VyBz 
0 


Finally the source vector S is given by 
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0 

0 

0 

-pg 

-pvzg 

0 

0 

0 


2.2 Characteristic Structure 

The numerical methods used to solve the MHD equations 
require the characteristic information from the Jacobian ma¬ 
trices associated to the system of equations § 


^.(Q)-gQ, 

As we can see, G and H can be obtained from an adequate 
index permutation of F. Therefore, Ay and Az have simi¬ 
lar characteristic structure as Ax. The set of eigenvectors 
for the Jacobian matrices has been well stablished and ex¬ 


tensively studied by several authors (Roe 1986 Brio & Wu 
Powell 1994 Ryu Jones 1995[ Balsara 2001). The 


1988 


corresponding eigenvalues for Ax, following the 8x8 eigen- 


system of Powell (Powell 1994 Powell et ah 1999) are given 
by 


^^1,8 — '^x i Cf, A2,7 — i Ca, ^3,6 — '^x i Cg, A4^5 — Vx, 

where Ca is the Alfven speed 



Cf and Cs are the fast and slow magnetoacoustic speeds re¬ 
spectively 



and a = . / — is the sound 

V ^ 

are found for Ay and Az. 


speed. Finally, similar expressions 


2.3 Numerical Methods 


We solve numerically the magnetized Euler equations given 
by the system 0 -§ on a single uniform cell centered grid, 
using the method of lines with a third order total varia¬ 
tion diminishing Runge-Kutta time integrator (RK3) p7i| 


& Osher 1989). In order to use the method of lines, the right 


hand sides of MHD equations are discretized using a finite 
volume approximation with High Resolution Shock Captur¬ 
ing methods. For this, we first reconstruct the variables at 
cell interfaces using the Minmod and MC linear piecewise 


second-order reconstructors, and the fifth-order weighted 


Essentially Non Oscillatory (WEN05) method (Titarev & 


Toro 2004 Harten A. et al. 1997| [Radice &: Rezzolla 2012 ) 
that uses a polynomial interpolation with a fifth order pre¬ 
cision for smooth distributions and third order at strong 
shocks. Even though the accuracy of the RK3 dominates 
over a high order accuracy reconstructor during the time in¬ 
tegration, we include WEN05 because it introduces a small 
amount of dissipation. 

On the other hand, the numerical fluxes are built with 
the Harten-Lax-van-Leer-Einfeldt (HLLE) approximate Rie- 
mann solver formula, using the left and right states of the 


primitive and conservative reconstructed variables (Harten 


P. et al. 1983 Einfeldt 1988). The HLLE approximate solver 


uses a two-wave approximation to compute the fluxes across 
the discontinuity at the cell interfaces, including those of the 
magnetic and hydrodynamic waves. Such waves move with 
the highest and lowest velocities that correspond to eigen¬ 
values of the Jacobian matrix. These numerical methods 
were implemented successfully in relativistic MHD, where 
we have studied different phenomena with fixed and dy¬ 
namic background space-times, including the presence of 


black holes (Lora-Clavijo 2014 Lora-Clavijo et al. 2013 


Lora-Clavijo et al. 2013b| |Lora-Clavijo Guzman 20T3 

Cruz-Osorio et al. 2012| [Lora-Clavijo Cruz-Osorio 2015| ) 


Additionally, in order to preserve the constraint ^ and 
avoid the violation of Maxwell’s equations (i.e. the devel¬ 
opment of unphysical results like the presence of magnetic 
charges) we use the flux constraint transport (flux-CT) al¬ 
gorithm (Evans & Hawley 1988 Balsara 200 Ij ). The flux-CT 
method is adapted to the conservative scheme computing the 
fluxes at cell vertex as an average of the fluxes at intercell 
boundaries and resulting in a cell-faced evolution equation 
of the components of the magnetic field. This scheme main¬ 
tains the divergence of the magnetic field computed at cell 
corners to round-off accuracy along the numerical evolution. 
The details of the numerical implementation can be found 


in (Lora-Clavijo et al. 2015) and other divergence constraint 


control methods are summarized in (Toth 2000). 


3 NUMERICAL TESTS 

We include a set of ID and 2D tests that involve the evo¬ 
lution with the gravitational field switched off. The solar 
tests include 2.5D and 3D tests involving the propagation of 
MHD waves in the solar atmosphere submitted to a constant 
gravitational field along the vertical direction. 


3.1 ID tests 


In order to illustrate how our implementation handles the 
evolution of a magnetized gas, in this subsection we present 
the standard shock tube tests in one direction. The ID tests 
are Riemann problems under various initial conditions and 
we compare the numerical results with the exact solution 
calculated in ( Takahashi Yamada 2Q13| Takahashi & Ya- 
mada 2014). 


We calculate the numerical solution of these tests us¬ 
ing the HLLE formula and the various limiters mentioned, 
in the domain [—0.5, 0.5] with N = 1600 identical cells, a 
constant Courant factor CFL = 0.25 for the first test, and 
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Test type 


7 


Brio — Wu 

Left state 5/3 

Right state 


1.0 

1.0 


1.0 

1.0 


0.0 

0.0 


0.0 0.0 

0.0 0.0 


0.75 

0.75 


1.0 

- 1.0 


Ryu&Jones lb 

Left state 
Right state 


Ryu&Jones 2b 

Left state 
Right state 


5/3 


1.0 

0.1 


1.0 

10.0 


0.0 

0.0 


0.0 

2.0 


Ryu&Jones 3b 

Left state 
Right state 


5/3 


1.0 

1.0 


1.0 

1.0 


- 1.0 

1.0 


0.0 

0.0 


0.0 

0.0 


0.0 

0.0 


1.0 

1.0 


0.0 

0.0 


5/3 1.0 1.0 0.0 0.0 0.0 3/V4^ 5/^4^ 0.0 

0.1 10.0 0.0 0.0 0.0 3 / V 47 r 2 /^/^ 0.0 


0.0 3/\/4^ 6/\/4^ 0.0 

1.0 3/V4^ l/Vi^ 0.0 


0.0 

0.0 


Table 1. Parameters for the various MHD ID Riemann problems. 




Figure 1. Test 1: Brio-Wu shock tube problem at t = 0.1. We 
show the numerical solution for p, Vx and By (points) and the 
exact solution (line). In this case we use the Minmod limiter. 


0.05 for the others. The reason to choose this value is that in 
order to compare among reconstructors used here, we use a 
value that can be handled by all of them. The initial discon¬ 
tinuity in this case is located at x = 0. In order to estimate 
the error between the numerical solution and the exact solu¬ 
tion, we calculate the numerical solution using four spatial 
resolutions Axi = 1/200, Ax 2 = 1/400, Axs = 1/800 and 
Ax 4 = 1/1600. For these tests we are using the 3D code 
with four cells along the transverse directions. The various 
parameters of ID tests are summarized in Table and we 
practiced this test for the three limiters mentioned. 


3.1.1 Test 1: Brio-Wu Shock Tube Problem 


The Brio-Wu MHD shock tube problem (Brio & Wu 1988) 


tests the capability of a code to handle waves of various 
types, shocks, rarefactions, contact discontinuities and the 
compound structures of MHD. A snapshot of the wave struc¬ 
ture at t = 0.1 of the solution is shown in Fig. The left- 
moving fast rarefaction wave from x ~ —0.2tox ~ —0.1 and 
the slow compound wave at x —0.03 are well captured. 
The contact discontinuity is captured at x 0.05. Finally 
the slow shock is captured at x ~ 0.15 and fast rarefaction 
from X ^ 0.32 to x 0.37. 


Figure 2. Test 2: Ryu & Jones lb shock tube problem at t = 0.03. 
We show the numerical solution for p, Vx and By (points) com¬ 
pared with the exact solution (line). Again we use the Minmod 
limiter. 



Figure 3. Test 3: Ryu & Jones 3b shock tube problem at t = 0.1. 
We show the numerical solution for p, Vx and magnetic field By 
(points) compared with the exact solution (line). 


3.1.2 Test 2: Ryu & Jones lb 


This test was introduced by (Ryu & Jones 1995). It is im¬ 


portant to check on the ability of the numerical scheme to 
handle a fast shock, a rarefaction shock, a slow shock, a 
slow rarefaction, and a contact discontinuity. A snapshot of 
the evolution is displayed in Fig. Our results show that 
fast shock at x ^ —0.11 and the slow shock at x ^ —0.075 
are well captured. In addition the contact discontinuity ap¬ 
pears at X —0.061 and the slow rarefaction is located 
dX X ^ —0.040. Finally the fast rarefaction is captured at 
X 0.34. 


3.1.3 Test 3: Ryu & Jones 3b 

This is a test containing magnetosonic rarefactions waves. 
The numerical solution shows only two strong and identical 
magnetosonic rarefr act ions. Our code was able to capture 
this structure and the results are shown in Fig. 
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3 . 1.4 Error estimates for the ID MHD tests 

In order to compare the numerical solutions using different 
limiters implemented in our code, we have calculated the 
error for each ID Riemann test described above, and the 
results are summarized in Table In all the tests, with the 
three reconstructors used, our runs lie within the consis¬ 
tency regime, that is, the error decreases when resolution is 
increased. In order to go further we check for the order of 
convergence. For this we calculate the Li norm of the error 
as compared with the exact solution. Since we always con¬ 
sider resolution factors of two, the order of convergence is 
given by log(errori/errori_i)/log(2), where i is the error 
calculated with the resolution Axi. Thus we carried out the 
tests with four spatial resolutions. For all the tests we ob¬ 
tained near first order convergence as expected, considering 
all the initial data start with discontinuities and in Table 
l^we show how the different reconstructors perform in the 
different tests. 


3.2 2D tests 

3.2.1 Orszag-Tang Vortex 


The Orszag Tang vortex system was first studied by (Orszag 


& Tang 1998). This is a well known problem for testing the 


transition to a supersonic 2D MHD turbulence, and to mea¬ 
sure the robustness of the code at handling the formation 
of MHD shocks. The important dynamics developed during 
the evolution is a test for the controller of the magnetic field 
divergence free constraint V ■ B — 0. 

The initial data are v = r'o(—sin(27ry), sin(27rx), 0), 
B = Ho(—sin(27ry), sin(47rx), 0), po = 25/(367r) and po = 
5/(127r), where vq — 1 and Bq — l/\/4^) with 7 = 5/3, are 
defined in such a way that the system develops a vortex. 
The initial data for p and p are calculated using averaged /3 
and the Mach number M (Dai & Woodward 1998), which 
are defined as 


/3 

M 


PO o 

< 

Cs 7Po' 


In this case /3 = 10/3 and M = 1. We carry out the 
evolution on the domain [0,1] x [0,1] that we cover with 
512 X 512 cells, and impose periodic boundary conditions. 

The evolution of the density shows the development of 
turbulence and the interaction of the shock within the do¬ 
main as shown in the snapshots taken at times t = 0.5 and 
t = 1.0 in Fig. [ 4 } which can be co mpared with the results 
from ATHENA ^tone et al. 2008 ). Our code maintains the 
constraint on values |V • B\ 10“^^, which is within the 


round-off error of the machine accuracy. In order to show 
how the fluid makes the transition to a supersonic turbu¬ 
lence, we calculate the Mach number defined as M = \v\^/c^s^ 
where \v\^ — where ^ is the speed of sound. 

Snapshots taken at times t = 0.5 and t — 1 oi M are shown 
in Fig. where we can see that Mach number increases 
in regions where the density decreases. More quantitative 
results are shown in Fig. where the pressure with vari¬ 
ous resolutions is measured on the line y = 0.3125 at time 


1 1.5 2 2.5 3 3.5 4 



0 0.2 0.4 0.6 0.8 1 


1 1.5 2 2.5 3 3.5 4 



0 2e-13 4e-13 6e-13 8e-13 1e-12 



Figure 4. Density for the Orszag-Tang vortex at times t = 0.5 
(top) and t = 1.0 (middle). We show |V • H| at t = 1.0 in the xy 
plane (bottom). 


t = 0.5. The results are shown for the Minmod reconstruc¬ 


tor and are in agreement with those in Fig. 11 of (Londrillo 


& Del Zanna 2000). Moreover, the results can be compared 


with those using FLASH (Fryxell et al. 2000 


2009 Lee 2013), ATHENA (Stone et al. 2008) and PLUTO 


Lee & Deane 


(Mignone et al. 2007). 
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Resolution 

MM 

MC 

WEN05 

MM 

MC 

WEN05 


Error 


Order of 

convergence 


Brio — Wu 

Axi 

O.lle-1 

O.lle-1 

0.82e-2 




Ax2 

0.66e-2 

0.66e-2 

0.45e-2 

0.74 

0.74 

0.86 

Ax3 

0.41e-2 

0.41e-2 

0.28e-2 

0.69 

0.69 

0.68 

Ax4 

0.25e-2 

0.25e-2 

0.19e-2 

0.71 

0.71 

0.55 


Ryu&Jones lb 


Axi 

0.32e-l 

0.27e-l 

0.24e-l 




Ax2 

0.22e-l 

0.15e-l 

0.13e-l 

0.71 

0.84 

0.88 

Ax3 

0.14e-l 

0.80e-2 

0.74e-2 

0.65 

0.90 

0.81 

Ax4 

0.86e-2 

0.50e-2 

0.45e-2 

0.70 

0.67 

0.71 


Ryu&Jones 3b 


Axi 

0.71e-2 

0.46e-2 

0.43e-2 




Ax2 

0.37e-2 

0.23e-2 

0.22e-2 

0.94 

1.0 

0.96 

Ax3 

0.19e-2 

0.12e-2 

0.12e-2 

0.96 

0.93 

0.87 

Ax4 

O.lOe-2 

0.67e-3 

0.79e-3 

0.92 

0.84 

0.60 


Table 2. Li norm of the error in density calculated with four resolutions and three numerical reconstructors. We also show the order of 
convergence between the different pairs of resolutions. 


3.2.2 The circularly polarized Alfven Wave Test 


This test problem was presented in (Toth 2000). The circu¬ 
larly polarized wave is an exact nonlinear solution to the 
MHD equations, which allows one to test a code in the 
nonlinear regime. This test is useful to measure the abil¬ 
ity of the numerical scheme to describe the propagation of 
Alfven waves. The initial data of this problem are p = 1, 
p = 0.1, v± = 0.1 sin(27rx||), i’ll = 0, B± = 0.1 sin(27rx||), 
Vz = 0.1 cos(27rX||), Bz = 0.1 cos(27rX||), where 
X|| = (x cos a+ y sin a) and a = 30° is the angle at which the 
wave propagates with respect to the grid. Here v± and B± 
are the components of velocity and magnetic field perpendic¬ 
ular to the wavevector. The Alfven speed is \va\ = B\\/p = 1, 
which implies that at time t = 1 the flow is expected to re¬ 
turn to its initial state. 


The domain [0,1/cos a] x [0,1/sin a] is covered with 
256 X 256 cells. In order to see that our code is capable 
to describe the propagation of circularly polarized Alfven 
waves along the diagonal grid, in Fig. (top) we show a 
snapshot of B± as function of i’ll taken at time t = 5 for 6 
different resolutions in the case of a travelling wave (r'y = 0) 
(top), in this case the diffusion of the wave is less for a higher 
resolution as expected. We also show a snapshot of Bz taken 
at time t = 5, where one can notice that the wavefronts 
remain planar despite the propagation occurs along a non¬ 
preferred grid direction. Again the code maintains | V ■ B\ 
10“^^ as shown in Fig.|^ In this test we use Minmod and the 


results can be compared with those obtained with ATHENA 
(Stone et al. 2008). 


3.2.3 Current Sheet Test problem 

This two dimensional problem was presented in 
) and is important because it shows 
to deal with magnetic reconnection driven solely 
by numerical resistivity. Additionally, the fact that magnetic 
reconnection produces high pressure gradients this test also 
serves to monitor the behavior of the code when dealing with 
these sharp gradients. 

The computational domain [—0.5,0.5] x [—0.5,0.5] is 
covered with 600 x 600 cells. We use the Minmod reconstruc¬ 
tor and impose periodic boundary conditions in all sides. We 
initialize the two current sheets as follows: 

^ ^ f Ho for |x| > 0.25 
^ \ —Ho otherwise^ 

where Ho = 1. The other magnetic field components H^, Bz 
are set to zero. The velocity is v = (i^o sin(27ry), 0, 0) with 
uo = 0.1. The density p = 1 and the gas pressure p = /3/2, in 
this case we choose 13 = 0.1. The problem is initialized with 
a magnetic field along the y direction, that switches signs 
twice in the domain. This system is then perturbed with a 
sinusoidal velocity function, which initially creates nonlin¬ 
ear polarized Alfven waves. These waves quickly turn into 


Stone 1995 


of the code 


Hawley & 
the ability 
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0 0.2 0.4 0.6 0.8 1 



Figure 5. Mach number M for the Orszag-Tang vortex at times 
t = 0.5 (top) and t = 1.0 (bottom) in the xy plane. The supersonic 
regions are located in the dark colored zones. 



Figure 6. Horizontal slices of p along y = 0.3125 for the Orszag- 
Tang vortex using four spatial resolutions using Minmod. This 
plot shows the consistency of the solution. 


magnetosonic waves. Eventually, numerical reconnection oc¬ 
curs between the field lines as shown in Fig. [^(bottom- left), 
which causes magnetic energy to be turned into thermal en¬ 
ergy. This process can be seen in the plots of Fig. which 
show regions of high gas pressure and low magnetic pres¬ 
sure, these regions appear in the form of bubbles which 
eventually merge with each other. Moreover, even though 
reconnection takes place, the constraint is kept under con¬ 
trol I V'B| 10“^^. The results can be compared with those 


using FLASH (Fryxell et ah 2000 Lee & Deane 2009 Lee 


2013) and ATHENA (Stone et ah 2008). 



Figure 7. Circularly polarized Alfven wave test. We show the 
perpendicular component of magnetic field as a function of 
at time t = 5, using six resolutions (top), Bz at time t = 5 
(middle) and |V • H| at t = 5 (bottom) in the xy plane. 
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Figure 8. Current sheet test. We show gas pressure p, magnetic pressure B^/2, magnetic field lines and |V • B\ at time t = 1.5 in the 
xy plane. 


3.3 Solar Tests 

In order to show that our code can handle solar physics 
problems in the gravitationally stratified solar atmosphere, 
we study the transverse oscillations of Alfvenic pulses in so¬ 
lar coronal loops using a 2.5 D model, and in the case of a 
3D configuration, we study the propagation of impulsively 
generated magnetohydrodynamic-gravity waves and vortices 
in the gravitationally stratified non-isothermal solar atmo¬ 
sphere. 


3.3.1 Transverse oseillations in solar eoronal loops 


Following the idea developed in (Del Zanna et al. 2005), we 
study impulsively generated Alfven waves propagating in a 
magnetic arcade, taking into account a stratified solar at¬ 
mosphere with height from the photosphere-corona region 
in order to study transverse oscillations in coronal loops. 
These oscillations are considered due to flare events gener¬ 
ally interpreted as standing kink oscillations. 

In this problem, the model of the coronal arcade is 
obtained from a solution of the static 2D MHD equations 


(Priest 1982), and is given by: 


Bx = Bq cos{kx) ex.'p{—kz)^ 

Bz — —Bosm{kx)ex.p{—kz), (10) 

where Bo = 40 G is the photospheric held magnitude at the 
footpoints X = ztLl2 and k = tt/L. We choose in our set up 
to use L = 50Mm. This is a force and current free magnetic 
held, which implies that the pressure gradient must simply 
balance gravity in the vertical direction, that is, according 
to ([^ 


-Vpe + Peg = 0. (11) 

Using the equation of state for a fully ionized hydrogen 
plasma p = TkspT we obtain the thermal pressure: 


p(z) = p(zo) exp 



( 12 ) 


where p{zo) is the pressure at photospheric level. Then the 
mass density is calculated using the equation of state ^ 
and gives p{z) = problem the temperature 

model for the solar atmosphere is given by a smoothed step 
function 
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T{z) = -Tear 


1 + dtc + (1 — dtc) tanh 


z — Zt 


(13) 


where dtc = TphotiTcor and Tphot = 6000 K is the temper¬ 
ature at the photosphere, Tcor = 1-2 x 10® K is the higher 
coronal temperature. These two regions are separated at 
Zt — 2 Mm by a transition region of width Zw — 0.2 Mm. 


With this profile we integrate (12) to obtain the pressure 
and therefore the density. 

Following (Del Zanna et al. 2005), in order to simu¬ 


late impulsively generated Alfven waves, we add a purely 
transversal velocity pulse Vy at initial time t = 0 , located at 
{xq,Zq)\ 


where Bo = 23 G, which is a value close to a typical magnetic 


field magnitude in the quiet Sun (Vogler & Schiissler 2007). 


As a result of applying the static equilibrium and the 
magnetic current-free conditions to equation the pres¬ 
sure gradient is balanced by gravity through equation ( 11 ). 
In addition, considering the fact that hydrostatic equilib¬ 
rium will only hold along the 2 : direction, the equilibrium 
gas pressure takes the form: 


Pe{z) “ 


[-f 


dz' 

W)j' 


(16) 


and the mass density 


_ Avq 

~ l + (r/ro)4’ 


(14) 


where r = {x — xo)^ + (z — with xo = 0 Mm and 
zq = 50 Mm, A = 0.1 is the normalized amplitude, vo = 1 
Mm/s is the typical Alfven speed in the corona and ro = 1 
Mm. Therefore the initial conditions for this problem are 
given hy : p = p(z), p — p(^), the velocity field v = ( 0 , 0 ), 

and B = (5^,0, B^). 

We solve the equations using the Minmod slope 

limiter, a Courant factor CFL—t).2 and 7 = 5/3. We set 
the simulation domain to [-25,25] x [0,1] x [0,50] Mm®, cov¬ 
ered with 400x4x400 cells, and we impose open boundary 
conditions in all faces. 

We consider the case of a symmetric pulse given by the 
equation ( |14[ ). In Fig.we show snapshots of the transverse 
velocity Vy at times t = 10 s, t = 40 s, t = 100 s and 
t = 200 s. At time t = 10 s we can see the propagation of two 
Alfvenic pulses with the same amplitude move downward 
following the fieldlines of coronal loop. At time t = 40 s the 
pulses have passed transition region where their amplitude 
decrease. The pulses move in regions with low density at 
time t = 100 s, so the pulses assume a elongated shape. At 
time t = 200 s the pulses hit the transition region twice and 
came back to the center of domain. In this time it is clear 
the spreading of the pulse, that goes together with a high 
amplitude decrease produce by reflection at the transition 
region. Even though the complex dynamics of the Alfvenic 
pulses, the V-B is maintained in the order of 10~^^ Tesla/km 
as is shown in Fig.(bottom). 


3.3.2 Propagation of magnetohydrodynamic-gravity waves 
in the 3D solar atmosphere 


Following (Murawski et al. 2013), we simulate the evolution 


magnetohydrodynamic-gravity waves and the vortex forma¬ 
tion in the magnetized three-dimensional (3D) solar atmo¬ 
sphere. Initially we assume the solar atmosphere to be in 


static equilibrium (jMurawski & Musielak 2010| 

Murawski & 

Zaqarashvili 2010||Murawski et al. 2011||Wolosz 

iiewicz et al. 

2014| IGhmielewski et al. 2014a| Ghmielewski et al. 2014b| 

Gonzalez-Aviles & Guzman 2015), and the initial equilib- 


rium magnetic field Be constant along the vertical direction 


Be = Boz, 


(15) 


Pe{z) = 


Pe{z) 

qHzY 


(17) 


where 


A(z) = 


ksTeiz) 

mg 


(18) 


is the pressure scale-height, and pref represents the gas pres¬ 
sure at a given reference level, that we choose to be located 
at Zref = 10 Mm, because in the temperature model above, 
this is the location at which the corona region starts. In this 
case we use p =1.24 in (1^ only to reproduce the results in 


(Murawski et al. 2013). 


The temperature field Te{z) is assumed to obey the 
C7 model ( Avrett fc Loeser 2008| , which is a semiempirical 
model of the chromosphere. This C7 temperature profile was 
interpolated into the 3D numerical grid, then we numerically 


integrated equation (16) in order to obtain Pe{z) and finally 
with this profile we use equation (17) to obtain pe{z). We 
show the temperature and density profile in Fig. |10| where 
the important density gradients can be observed. In addi¬ 
tion we show the plasma /3 to see how the magnetic pressure 
dominates over gas pressure in the solar atmosphere. 

In order to simulate the propagation of 
magnetohydrodynamic-gravity waves, we perturb the 
static solar atmosphere with a Gaussian pulse, either in 
the horizontal component of velocity, Vx, or in the vertical 


component Vz. Following (Murawski et al. 2013), the initial 
conditions are given by 


_ x‘^+y‘^ + (z-ZQp 

Vx{x,y,z,t = 0) = Aye , (19) 

_ x‘^+y‘^ + {z-ZQp 

Vz(x,y,z,t = 0) = Aye , (20) 

where Ay is the amplitude of the pulse, zq is the inital po¬ 
sition over the vertical direction and w denotes the width 
of the pulse. In this case we use the values Ay — kms“^, 
It; = 100 km and zo = 500 km. Therefore the initial condi¬ 
tions for this problem are given by: p = Pe{z), p = peiz), 
the velocity field v = {vx,0,Vz) depends on the direction of 
perturbation, and B = (0,0, Bo). 

We solve the equations 0-0 using the HLLE Rie- 

mann solver, the Minmod slope limiter, a Gourant fac¬ 

tor CBL=0.2 and 7 = 1.4. We set the simulation do¬ 
main to [-0.75,0.75] X [0,3] X [-0.75,0.75] Mm®, covered with 
64x128x64 cells for the case of horizontal perturbation and 
100x200x100 cells for vertical perturbation. At the top and 
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Figure 9. Transverse velocity ^^/(km/s) at times t = 10 s, t = 40 s, t = 100 s, t = 200 s and V • B (Tesla/km) at t = 200 s. 


bottom faces of the numerical domain we impose fixed in 
time boundary conditions, whereas at the remaining four 
sides we apply open boundary conditions. In our simula¬ 
tions we choose the scale factors in cgs units as specified in 
Table These units are obtained choosing a length-scale 
/o, plasma density scale po and magnetic field scale Bo, this 
fixes the unit of time in terms of Alfven speed as follows: 


Bo — r’oy/MoPo, (21) 

where po = 47r x 10“^ N/A^ is the magnetic permeability 


and po = 10 kg/m^ is the plasma scale at the level of the 
solar corona. 


3.3.3 Vertical perturbation 

Vertical perturbations in the sun could be generated by p- 
modes (Jain R. et al. 2011), which are pressure or acoustic 
modes with frequencies greater than 1 mHz. This kind of 
perturbation is given by the component Vz in equation (20). 
The initial pulse triggers longitudinal magnetoacoustic- 
gravity waves. At the initial point where the wave is 
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z(Mm) 


Figure 10. Temperature, density and /3 profile for the C7 model. 


Table 3. Units used in this paper 


Quantity 

cgs units 

lo 

10® cm 

VO 

10® cm-s“^ 

to 

1 s 

PQ 

10-15 gr.cm-® 

Bo 

11.21 G 


launched, the gas pressure is greater than the magnetic pres¬ 
sure, which implies that /3 > 1. This excites fast and slow 
magnetoacoustic-gravity waves which are coupled. The ini¬ 
tial pulse splits into two propagating waves at time t — 37.5 
s shown in Fig. (top). At time t = 100 s, the two mag- 
netocoustic waves move upwards until z 1.4 Mm, in 
the same way as in Fig. 3 of (Murawski et al. 2013). The 
waves continue moving until they reach the transition re- 


z 





Figure 13. Temporal snapshots of streamlines at t = 37.5 s (top) 
and t = 150 s (bottom) for the case of vertical perturbation. Color 
represents the vorticity magnitude. 


gion. This vertical perturbation produces the propagation 
of magnetoacoustic-gravity waves in a symmetric way, as 
shown in Fig.[^at the xz and yz planes. In these two planes 
the propagation is seen in the same way. We also show the 
propagation as seen from above, on the xy plane in Fig, 
Even though the structure of the magnetoacoustic-gravity 
waves in the solar atmosphere looks complex, the value of 
|V • at time t = 150 s in the xy plane remains low as is 
shown in Fig. |12[ 

We show snapshots of the fluid streamlines at times 
t = 37.5 s and t = 150 in Fig. These streamlines reveal 
vortices specially at time t = 37.5 s, which are produced by 
the vertical perturbation Vz and are located at 1 Mm over 
the photosphere. 


3.3.4 Horizontal perturbation 


Horizontal perturbations in the solar atmosphere can be gen¬ 
erated due to granular cells as a result of convection to¬ 
ward the photosphere. In this case we simulate these waves 


by exciting the Vx component in equation (19). This ini¬ 
tial pulse excites all MHD waves which are affected by the 
gravitational constant field. Since early times we can see the 
magnetoacoustic-gravity waves propagating upwards in the 


xz plane as ligth yellow wave-fronts in Fig. 14 (top). At time 
t — 150 s we can see the magnetoacoustic-gravity waves light 
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t = 37.5 s 


t = 100 s 


t = 150s 


Vz{km/s) in yz plane 



Vz{km/s) in xz plane 
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Figure 11. Snapshots of the vertical component of velocity Vz{km/s) at times t = 37.5,100,150 s in the yz and xz planes. 
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t = 37.5 s 
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Figure 12. Snapshots of the vertical component of velocity Vz{km/s) at times t = 37.5,100,150 s and |V • S| at t = 150 s in the xy 
plane. 


yellow patches reaching the level z = 1.75 Mm as shown in 


the top of Fig. 6 of (Murawski et al. 2013). At time t = 250 


s, the waves reach the transition region, which causes the 
appareance of a light blue patch in Fig. but later on, 
they are reflected back. In order to see the propagation of 
the waves from different planes, we show Vx projected on the 
xz and yz planes in Fig. (bottom). Unlike the previous 
test, in this case the velocity is not symmetrical. Moreover, 
the propagation of the waves seen from above show how 
the wave-fronts change their amplitudes until the reflection 
happens, this is shown in the top panel of Fig. |15| As in the 
previous case we show a snapshot of | V • 5| at time t = 150 


s in Fig. 15 (bottom). 


In this case we also show snapshots of the fluid stream¬ 
lines at times t = 50 s and t = 150 s in Fig. In the same 
way as in the case of the vertical perturbation, the hori¬ 
zontal perturbation generates vortices, but in this case the 
vorticity is more complex. It has been proposed that vortices 


are triggered by convective or granular motions (Murawski 


et al. 2013). 


4 FINAL COMMENTS 

We have presented a new code that solves the ideal MHD 
equations intended to model phenomena in the solar atmo¬ 
sphere. 

We have included standard ID Riemann problem tests 
and shown that the order of convergence lies within the 
expected range corresponding to the methods we use. We 
conclude from our analysis, that MC and WEN05 perform 
better than the Minmod limiter. 

We have also presented standards 2D tests including 
the Orszag-Tang vortex and a current sheet test problem 
showing magnetic reconnection. In these tests we have kept 
the magnetic field divergence free constraint under control, 
in all the tests near round-off error. We notice that, accord¬ 
ing to previous experience, due to the dissipation, Minmod 
performs well and the MC and WEN05 limiters allow more 
refined structures in general. 

The solar tests have been oriented to study wave prop¬ 
agation in a gravitationally stratified solar atmosphere. We 
track the transverse oscillations of Alfvenic pulses in solar 
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t = 50 s 


t = 150 s 


t = 250 s 


Vx{km/s) in xz plane 
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Figure 14. Snapshots of the horizontal component of velocity Vx{km/s) at times t = 50,150, 250 s in the xz and yz planes. 
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t = 50 s 
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Figure 15. Snapshots of the horizontal component of velocity Vx{km/s) at times t = 50,150,250 s and |V • S| at t = 250 s in the xy 
plane. 


coronal loops in a 2.5D model, in adittion we study the prop¬ 
agation of MHD-gravity waves and vortices on a 3D solar at¬ 
mosphere. These results were compared to previous results 
in the literature. We have also shown that the constraint has 
been kept under control. 

We have decided to use unigrid mode instead of mesh 
refinement in order to avoid the use of artificial viscosity at 
boundaries between refinement levels, which is used to damp 
out spurious oscillations. We also want to avoid the carbun¬ 
cle effect at corners of refined grids. In order to carry out 
production runs with high resolution, we have parallelized 
our code. 

Finally, the code is public and available under request 
at http://www.ifm.umich.mx/~cafe 
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